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Abstract 

Cellular neural circuit and networks consisting of interconnected neurons and glia are ulti¬ 
mately responsible for the information processing associated with information processing in the 
brain. While there are major efforts aimed at mapping the structural and (electro)physiological 
connectivity of brain networks, such as the White House BRAIN Initiative aimed at the devel¬ 
opment of neurotechnologies capable of high density neural recordings, theoretical and compu¬ 
tational methods for analyzing and making sense of all this data seem to be further behind. 
Here, we propose and provide a summary of an approach for calculating effective connectivity 
from experimental observations of neuronal network activity. The proposed method operates 
on network-level data, makes use of all relevant prior knowledge, such as dynamical models of 
individual cells in the network and the physical structural connectivity of the network, and is 
broadly applicable to large classes of biological and non-biological networks. 
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1 Introduction 


Cellular neural circuit and networks consisting of interconnected neurons and glia are ultimately re¬ 
sponsible for the information processing associated with information processing in the brain. While 
the entire network can range in size from a few hundred neurons, as in the case of the nematode 
worm C. Elegans , to several hundred billion in the human brain, distinct functions and subtasks 
are presumably carried out by the interaction within and between relevant microcircuits. Within 
a neuronal network, there are two organizational topologies or connectivity classes: structural and 
dynamic. Structural connectivity describes the physical locations and connections between cells 
in the network, while dynamic connectivity is transient and changing and describes how one cell 
affects another. It is a subset of the physical connectivity [2j. The neuroscience literature dis¬ 
tinguishes two types of dynamic connectivity, functional connectivity, which describes statistical 
but not casual correlations between cells (or populations of cells and even entire brain regions), 
and effective connectivity, which is stronger and assumes casual dynamic connectivity within the 
network H 12 Jl3] . In a neurobiological context, this means that cells that are physically connected 
need not necessarily signal each other. While there are major efforts aimed at mapping the struc¬ 
tural and (electro(physiological connectivity of brain networks, such as the White House BRAIN 
Initiative aimed at the development of neurotechnologies capable of high density neural recordings, 
theoretical and computational methods for analyzing and making sense of all this data seem to 
be further behind. And while the experimental and technical challenges of imaging and mapping 
structural connectivity are tremendous, the objectives and goals are clear, and continued technical 
advances have ensured continued progress. Theoretically though it is much less obvious how to 
analyze, model, and use this data for the purposes of understanding brain function. One logical 
first step is to ask what is the effective connectivity associated with a specific imaged or recorded 
pattern of dynamic activity in a measured network or neural circuit with a given structural con¬ 
nectivity. Achieving this is at the forefront of theoretical and computational neuroscience, and 
certainly not a trivial task. Theoretically there are open mathematical problems that need to be 
solved in order to be able to accomplish this. Here, we propose and provide a summary of an 
approach for calculating effective connectivity from experimental observations of neuronal network 
activity. It builds on a neural signaling dynamic framework we have previously published [?]. The 
proposed method operates on network-level data, makes use of all relevant prior knowledge, such as 
dynamical models of individual cells in the network and the physical structural connectivity of the 
network, and is broadly applicable to large classes of both biological and non-biological networks. 

Several relatively recently developed methods have attempted to tackle the problem of func¬ 
tional connectivity estimation, approaching the problem from different perspectives. Eichler and 
colleagues developed a statistical-based approach that works with acyclic networks and computes 
direction and polarity of network connections |7|. Another method optimizes both neuronal and 
connectivity parameters of a simple linear integrate and fire deterministic model by comparing 
simulated output to experimental data, in this case being artificially generated spike trains 10 
. The method was validated for small networks of five neurons or less, and can handle feedback 
connections, though the estimated connection strengths are only reliable for indicating polarity and 
not relative synaptic strength. deFeo’s method uses three successive optimization steps to estimate 
a full state-space reconstruction of a measured signal, fitting of a local nonlinear dynamical model 
to the reconstructed signal, and then estimating a linear model to the interactions of the individual 
local models [5]. Finally, work by Eldawlatly, Zhou, Jin and Oweiss |8j, employs dynamic Bayesian 
networks for identification of connections of small (N = 10) networks using a Poisson spiking model. 
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More recently, Abarbanel and colleagues have described a powerful parameter estimation method 
optimized for tracking and predicting output voltages from unobservable internal state variables in 
models of dynamic neuronal signaling and coupled networks 00 However it remains unclear how 
their methods apply to and could take advantage of geometric networks where spatial information 
in addition to functional electrophysiological information is preserved. 

While these approaches are successful on the test data sets used for their validation, our proposed 
approach is broader and more applicable to experimentally derived observations. Our approach is 
model agnostic, meaning that it can work with many different neuronal and astrocytic individ¬ 
ual cell dynamics models. The only constraint is that the dynamics are observable, though this 
constraint is valid for any type of estimation. Our approach works with any type of network con¬ 
nectivity, meaning that a network can have recurrent connections (cyclic or acyclic). Finally, this 
approach takes advantage of the linear summation of signals into a cell to compute both functional 
connections and signaling delays which are present in biological networks, providing a novel ap¬ 
proach to estimating causal functional relationships in neuronal networks. We can also derive the 
minimum experimental data collection requirements, both in signal quantity and diversity, needed 
for mapping of a network of arbitrary size. 

Our previous work on network mapping methods and algorithms depends on stochastic non¬ 
linear identification methods that take advantage of high performance massive parallel graphics 
processing unit (GPU) computing architectures [?]. Although this approach is powerful, because 
of the dependency on stochastic numerical methods it can show significant variation in accuracy 
as a function of the non-linear optimizer used or the available computing power. By accuracy, we 
mean how close a cost function can be minimized in order to recover the real underlying effec¬ 
tive connectivity of the network for a given transient signaling event. While on-going research in 
both non-linear optimization and GPU computing will ensure that our current methods will only 
improve, a deterministic algebraic identification approach to appropriately constrained problems 
would ensure that as long as the constraints of any such method are satisfied, the accuracy of the 
mapped network will always be within arbitrary error bounds. 


2 Neural Cell Signaling in Geometric Networks 

We very briefly summarize our geometric network neural signaling framework. Uniquely, this frame¬ 
work takes advantage of the physical structure and connectivity of a neuronal network in addition 
the individual (electrophysiological) internal neuronal dynamics to determine and analyze the net¬ 
work dynamics. We refer the reader to [2] for complete details. We define a network as the graph 
G = ( V , E) to be a collection of vertices V and connectivity information E between vertices. The 
network is composed of N vertices, defined both geometrically (physical position) and dynamically 
(time variant state space): 

V = {xj,Yj(t): j = l,...,N} (1) 

Here x ? denotes the physical position in space of the vertex, while y j(t) is the state vector of the 
vertex at time t. Connectivity information is composed of a connective strength uj t j and time delays 
Tij , defined for all unique and directed pairs of vertices in the network: 

E = {uJij , Tij : i , j = 1,..., N, * 7 ^ J} (2) 

The connective strengths are scaled such that — 1 < Wy < 1 and delays are non-negative: 0 < Tij. 
We describe and validate the full framework in [?]. 
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Within a vertex, the time course of the state vectors is computed from the time map Hy between 
evenly spaced and sufficiently close time points, defined as follows: 


y j(t + A t) = Hj ( yj ( t ), sy (t), uj (t)) (3) 

Here, Uy (t) are the known external inputs to the system, and sy ( t ) represents the internal network 
signaling within the network. Relatively little restriction is placed on the functional form of Hy; 
it can be linear or nonlinear, stochastic or deterministic. Hy must be observable, meaning that all 
the states of yy can be estimated from observations of a subset of y r Several approaches exist for 
state estimation or filtering: particle filtering 16|, Kalman filters [9], and variational approaches 1 . 
The choice of estimation algorithm will depend on the individual cell models, and a full description 
of each is beyond the immediate scope of this paper. 

Intra-vertex network signaling must be of the form: 

N 

s j(t)= Y W/.'h'y-^ ~j)) ( 4 ) 

Under this assumption, the incoming signal to a vertex j is the sum of weighted functions of the 
delayed states of other vertices in the network. is the transmission function that operates on 
the state vector of another vertex. 

Finally, in similar fashion to the transmission function, we define an observation entity zy(t) as 
some function of the state vector y ? (f): 


z j(t) = F j(yj(t)) (5) 

z j(t) describes what is experimentally observable about the vertex j at time t. The state transition 
and observation functions Hy and Fj , must be observable so that estimation of signaling quantities 
is possible. 


Estimation with known or no delays 

When the delays are known or non-existent, equation [4] can be written as: 

N 

s j{t) = Y, (6) 

*=M#i 


where ry(t) = gi (y»(t— t^)). Under the system observability assumption, we can estimate both 
s j(t) and r’i(i) from the observed activity of individual cells z j(t) at each of the T time points, the 
functional connective weights oj,y can be obtained using some simple algebraic manipulation. We 
can expand the above equation to matrix notation as follows: 


Sj( 0) 

Sj(A t) 


n(0) r 2 (0) 

ri(At) r 2 (At) 

?*jv(0) 

r N (At) 


Uij 

UJ 2 j 

_ Sj((T — l)At) _ 


_n((r-l)Ai) r 2 ((T-l)At) • 

■ r N ((T — l)At) 


U>Nj 


4 








The above equation can equivalently be written as: 

Sj = R jUij (8) 

The vector Sj has dimension Txl, the matrix Rj is T x (N — 1) and the vector ujj, containing all 
the incoming functional weights into vertex j is thus (N — 1) x 1. For a unique solution to uij, two 
criteria must be met: T > N — 1 (ideally T N — 1), and Rj must have full rank (rank N — 1), 
meaning that no two rows of Rj can be linear multiples of each other. Since equation [8] is linear, 
any standard method of factoring out ojj can be used. For example, one can multiply both sides 
by the pseudo-inverse of Rj, isolating cjj to the right hand side. 


Estimation with unknown delays 


The more interesting and difficult case is when the delays of equation [4] are unknown. To address 
this problem, we first take the Fourier transform of [3] as follows: 


N 


hi 

CO 

II 

hi 

; Y W/f 
[i=M*7 J 

(9a) 

N 

s j( s )= Y 

co ljri (s)e- 2 ^ a 

(9b) 


In the frequency domain, equation |9b| employs the shift theorem, where a shift in the time domain 
of Tij results in a multiplication in the frequency domain by e “ 27rT ;j s . This property allows us to 
effectively remove the unknown delay Tij from the argument of the r* function and express it as the 
argument of the natural exponent. Since the exponent in |9b| has only an imaginary argument, the 
equation can be expressed as follows: 


N 

s j( s ) = Y^ Uijn(s) (cos(—27TTys) + isin(—27rr^s)) 

We can now separate the real and imaginary components of the above, as 

N 

Re[Sj-(s)] + ilm[sj(s)] = ^ (Re[r*(s)] + «Im[rj(s)]) (cos(—27rr^-s) + «sin(— 2nTijs)) 


( 10 ) 


( 11 ) 


The above complex-valued equation can be separated into the real and imaginary parts, into two 
real-valued sets of equations: 


N 

Re[sj(s)] = ^ ujij (Re[r - i(s)] cos(— 27 TTijs) — Im[rj(s)] sin(— 2nTijs)) 

N 

Im[sj(s)] = ^ Wjj (Re[rj(s)] sin(— 27 TTjjs) + Im[rj(s)] cos(— 27 TTjjs)) 


(12a) 


(12b) 
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We now have have a system of real-valued equations with and as unknowns and known 
values of s j(s) and rt(s) evaluated at every frequency s. The system is non-linear but can be solved 
using the Levenbergh-Marquardt solver (see 11 for description). Briefly, one can rewrite equations 


12 as a series of functions that equal 0, and find that values of 0 Jj : j and Tjj that make the above 


equations zero or bring them closest to zero. 

As an example, Fig. |T] shows a simulated example of transmitted signals in a network, and the 
measured summed and delayed activity of one vertex. Once the transmitted signals are estimated, 
the individual contribution (i.e. effective connectivity) and delays from each transmitted signal 
r,(f) into the compound signal s j(t) can be estimated by solving equation in 12 The results are 
shown in Fig. [2] In this example, we estimate 20 incoming functional connections and transmission 
delays into a vertex. To generate the synthetic data, a random spike train was generated for each 
cell as the transmitted signal, as shown in figure [Tj4 . Dynamic connections were randomly selected 
in a uniform distribution between -0.5 and 0.5, and delays were randomly chosen between 0 and 
400 time steps. The incoming signals were weight-modulated and delayed to produce the incoming 
signal Sj(t) shown in the bottom panel of Fig. 1. Using the transmitted signals and the incoming 
signals, our algorithm was able to fully recover the original weights and delays that were used to 
generate the synthetic data (without the algorithm being aware or told what those were, of course). 


Discussion 

This approach rests on the fundamental assumption that signals into a cell are summed and sub¬ 
sequently drive the dynamics of the cell in the network. This assumption generally holds true for 
most neuronal and astrocytic networks, where post-synaptic currents and extra-cellular ATP are 
the signaling quantities that sum up and drive the cell dynamical systems. From this assumption, 
we were able to express signaling in matrix form and perform the inversion with trivial matrix 
factorization (equation |8j) . What is apparent from this equation is that the number of time points 
collected from each cell has to be greater than the number of cells in the network for a unique 
solution to be computable. Additionally, the outgoing signals of other cells in the network must be 
different so that no two cells have the same exact dynamics, so that no two rows of the matrix Rj 
are identical, making Rj singular and non-invertible. When delays are present and significant, the 
problem is tractable, though requires a non-linear least squares solver to compute. An additional 
data requirement is that the individual outgoing signal differ in the frequency domain as well as in 
the time domain, meaning that two signals cannot be time-shifted copies of each other for a unique 
delay to be computable. Overall, these are necessary conditions, with sufficiency only possible for 
the linear, no-delay case. In the non-linear, delays case, only necessary conditions can be met; the 
non-linearities in the summation of sines and cosines make the problem more difficult, mandating 
use of an optimizer. Indeed few non-linear root-solving problems have sufficiency criteria for finding 
globally optimal solutions, as in the linear case. While not unique to this problem, nonlinear root 
finding is the major theoretical obstacle in this algorithm. 

Experimentally, the estimation approach we introduce here is tailored toward a mixed voltage 
and calcium recording experimental setup for neuronal functional connectivity. Membrane voltage 
recordings from a few neurons at a time can be used to estimate the applied synaptic current, while 
calcium fluorescence observations of the complete network provides estimates of the spike times and 
thus the possible transmitted signals. Combined, the incoming functional connections and delays 
are computed for the neurons whose voltage time courses are being recorded. 
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Figure 1: Incoming and Compound Measured Signals. Measured transmitted (top panel) 
and compound signals into a vertex (bottom panel). The top figure shows the network transmitted 
signaling activity, while the bottom shows the driving compound signal into one vertex. Both of 
these are estimated from the observed activity of the individual vertices. 


7 




















































































Functional Connectivity Estimation Transmission Delay Estimation 




Figure 2: Connectivity and Delay Estimation. Connectivity and Delays estimation for signals 
in Fig. 1. The system was solved using the Levenberg- Marquardt algorithm for nonlinear systems 
of equations. In this example, we estimated exactly the connective strengths and delays of the 
signals into the observed compound signal. 
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